Magnetic spectrum of trigonally warped bilayer graphene — semiclassical analysis, 
zero modes, and topological winding numbers 
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We investigate the fine structure in the energy spectrum of bilayer graphene in the presence of 
various stacking defaults, such as a translational or rotational mismatch. This fine structure consists 
of four Dirac points that move away from their original positions as a consequence of the mismatch 
and eventually merge in various manners. The different types of merging are described in terms of 
topological invariants (winding numbers) that determine the Landau-level spectrum in the presence 
of a magnetic field as well as the degeneracy of the levels. The Landau- level spectrum is, within 
a wide parameter range, well described by a semiclassical treatment that makes use of topological 
winding numbers. However, the latter need to be redefined at zero energy in the high-magnetic-field 
limit as well as in the vicinity of saddle points in the zero-field dispersion relation. 

PACS numbers: 73.43.Nq, 71.10.Pm, 73.20.Qt 



I. INTRODUCTION 

Graphene research has stimulated many fields of 
condensed-matter physics during the last years.— One of 
the most remarkable of these fields is certainly the topo- 
logical description of electronic energy bands, such as in 
the context of topological insulatorsi^ Indeed, the low- 
energy electronic properties of a single graphene layer 
are determined by two particular band-contact points at 
the corners K and K' of the first Brillouin zone, with 
a linear dispersion relation (the so-called Dirac points). 
These Dirac points are associated with a topological 
Berry phase that stems from the winding of the phase 
in the electronic wave function on closed paths around 
these points - the Berry phase is then tt times this wind- 
ing number i^ Prominent consequences of this Berry phase 
are the absence of backscattering in the case of long-range 
disorder—, Klein tunnelingf^i^ and a particular form of the 
Landau level (LL) spectrum in the presence of a magnetic 
field, with a topologically protected zero-energy levelj^ 

Bilayer graphene, that is obtained from an AB stacking 
of two graphene layers, has an even richer band structure, 
also from a topological point of view, than monolayer 
graphene. Most of its electronic properties have success- 
fully been described in the framework of two parabolic 
bands with opposite curvature that touch each other at 
the Fermi level. As compared to monolayer graphene, 
the winding number associated with these band-contact 
points is twice as large4 This gives rise to a two-fold 
orbital degeneracy of the zero-energy level in the pres- 
ence of a magnetic field, in addition to the four-fold 
spin-valley degeneracy, and thus to a particular series 
of Hall plateaus that have been observed in quantum- 
Hall measurements 1^ However, this picture is only ap- 
proximately valid in an intermediate energy range (above 
~ 10 meV), whereas subordinate hopping terms yield a 
fine structure in the energy spectrum (called "trigonal 
warping" ) , in the form of four Dirac points with lin- 
ear dispersion, at lower energiesi^ This transition from 
four Dirac points, with unit winding numbers, to the 



parabolic regime with a winding number of 2 may be 
viewed as a finite-energy Lifshitz transitior>i^ between 
disconnected Fermi pockets at low energies and a simply 
connected Fermi sea (per valley) at higher energiesJ^ In 
contrast to earlier experiments on bilayer graphene, to- 
day's availability of high-quality samples allows one to 
probe now this low-energy regime in which quantum-Hall 
measurements indicate the presence of additional Dirac 
points^ and it is noteworth to mention that indica- 
tions of Lifshitz transitions have been found in cyclotron- 
resonance measurements in graphite J^ 

The fine structure of the energy spectrum of bilayer 
graphene is also interesting from the point of view of 
stacking defaults, such as a displacement, strain or a 
twist with respect to perfect AB stackingiiirJ^ In this 
case, the low-energy dispersion is modified and two or 
more of the Dirac points may easily mergeJ^iii This 
needs to be contrasted to Dirac-point merging in mono- 
layer graphene that has been extensively studied on the 
theoretical leveU^"— but that is difficult to achieve ex- 
perimentally due to an enormous strain required)^ Fur- 
thermore, moderate stacking defaults may allow for the 
systematic study of merging transitions that fall into two 
distinct topological classesii - whereas the merging of 
Dirac points with opposite winding numbers yields a gap 
in the band structure, that of Dirac points with the same 
winding number maintains the band-contact points. This 
difference has direct consequences for the LL spectrum, 
namely the zero-energy level. Whereas in the former case 
of merging Dirac points with opposite winding number, 
the twofold degeneracy of the zero-energy level is lifted^^ 
it is topologically protected in the latter caseJ^ 

Here, we investigate the different merging transitions 
that one may encounter in bilayer graphene with a stack- 
ing default, within a continuum model that has been 
used both in the description of bilayer graphene with 
a mismatch described by a translation between the lay- 
ers or under strai n^^'^^ as well as in that of a twisted 
bilayer jii This continuum model, which goes beyond the 
linear Dirac-point approximation, may be viewed as a 



continuum model of the second generation)^ In addi- 
tion to the merging transition between Dirac points of 
opposite winding number, we discuss in detail the triple 
merging of three Dirac points that has been investigated 
in previous theoretical works>i^i^i21 This triple merging 
happens to be unstable in the sense that it only occurs 
in the framework of a displacement or strain in a high- 
symmetry axis of the lattice - a slight deviation from 
such an axis splits the triple- merging into a usual merg- 
ing transition of two Dirac points in a first step, followed 
by merging with the remaining Dirac point in a second 
step. As compared to previous studies of the LL spec- 
trum for trigonally warped bilayer graphene^^^ we pro- 
vide in the present paper a detailed semiclassical analysis 
of the spectrum. This analysis is based on winding num- 
bers that allow for a transparent understanding of LL 
degeneracies and zero modes. Furthermore, we investi- 
gate quantum corrections beyond the semiclassical limit. 
Both at zero energy in the high-field limit and in the 
vicinity of saddle points in the dispersion relation, these 
corrections are relevant because they blur the semiclas- 
sical trajectories and thus call for a modification of the 
description in terms of winding numbers. This allows for 
an understanding of the change in the LL degeneracy at 
zero energy and in the vicinity of saddle points in the 
dispersion relation. 

The paper is organized as follows. In Sec. |lTl we 
present the continuum model that accounts for the differ- 
ent stacking defaults in bilayer graphene and discuss the 
band structure and zero-field merging transitions in the 
several limits. Furthermore, we characterize the merg- 
ing transitions in terms of topological winding numbers. 
Section IIIII is devoted to the LL spectrum associated 
with the different merging transitions. The spectrum is 
obtained within a numerical solution of the quantum- 
mechanical eigenvalue equation (Sec. IIII Ap and ana- 
lyzed in the framework of a semiclassical treatment (Sec. 
IIIIBp . Topological aspects of the LL spectrum are dis- 
cussed in Sec. llVi and a detailed discussion of the dif- 
ferent aspects of the spectrum may be found in Sec. |Vl 
before we present our conclusions (Sec. IVip . 
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FIG. 1. (Color online) Lattice Structure of the Bilayer 
Graphene. (a) Atomic structure around an elementary cell 
for the Bernal configuration with {A,B) atoms (first layer) in 
blue and {A, B) atoms (second layer) in red. Hopping param- 
eters t, t± and i' are also pictured. Figures (b), (c) and (d) 
depict a planar lattice geometry for the Bernal, slided and 
twisted bilayer, respectively. 



A. Tight-Binding Approach 

For perfect AB-stackingf§. the tight-binding approxi- 
mation yields a four-band Hamiltonian that may be writ- 
ten in the {A, B, A, B) basis 



H(k) 



/ i7(k) i'7*(k) 

t7*(k) t_L 

t± t7(k) 

V t'7(k) t7*(k) 



(1) 



where 



II. BAND STRUCTURE OF DEFORMED 
BILAYER GRAPHENE 

Bilayer graphene harbors different stacking geometries, 
pictured in Fig. [1] among which the energetically most 
favorable is the Bernal (AB) configuration in which a 
B sublattice atom of the first layer sits on top of an A 
sublattice atom of the second layer [Fig. [l{a)]. This par- 
ticular ordering is naturally observed in graphite, as well 
as for synthesized bilayer graphene. Other stackings may 
be described in terms of a rotation defaull— and a dis- 
placement vector-i^^ and may be observed in graphene 
samples, such as for example in epitaxial graphene on 
the C-face of the SiC crystal.— We also consider strain 
constraints along both layers j^^j^S 



7(k) = -(l 



^ikai 



(2) 



and ai , a2 are elementary vectors of the triangular Bra- 
vais lattice 



ai,2 = -{±Viex + 3ej^) 



(3) 



and a = 0.142 nm is the distance between neighboring 
carbon atoms in the same layer. The hopping parameters 
can be experimentally evaluatedj^S and one obtains the 
hierarchy 

i(-3eV)>ti(-0.4eV) > t'(- 0.3 eV), (4) 

where t represents the hopping between p^ orbitals of 
nearest-neighbor carbon atoms within the same layer, t±^ 



the perpendicular hopping amphtude between a B sub- 
lattice atom of one layer and the A atom of the other 
layer, and t' is the transfer integral from an A site of 
one layer to the nearest B sites of the other layer [see 
Fig. [Ija)]. All other orbital overlap may be neglected for 
energies larger than 1 meV.— 



B. Low-Energy Hamiltonian 

In the small-wavc-vector limit (|q| <C 1/a), Hamilto- 
nian ([1]) may be expanded around a, K ot K' corner of 
the hexagonal Brillouin zone situated at the positions 
±K = ±4:Trex/3^/3a, modulo a reciprocal lattice vec- 
tor. One has then t7(±K + q) ~ vpi^Qx — iQy), in 
terms of the Fermi velocity vp = 3ta/2 {h — 1 hence- 
forth) and q -^ K. Furthermore, for energies lower than 
i_L, only two bands are relevant, and they may be de- 
scribed with the help of an effective two-band continuum 
Hamiltonian^ 



n 



K 



Here, \b\ 



^t2 

-K^ 



TT 
Tft 



— T~Lh + T~Lc 



(5) 



14: /mo, in terms of the bare elec- 



tron mass Too, c = vpt' /t sa 10^ m/s, and tt = q^ + iqy 
is the complex momentum operator in the continuum 
limit (that changes as tt ^- — tt^ when K — ^ K' = — K). 
The Jii, term in Eq. ([5]) is dominant for energies higher 
than ~ 10 meV and lower than t± ^ 0.4 eV. In the ab- 
sence of the term 'He, it enforces a quadratic dispersion 
around the band-contact points at K and K'. For en- 
ergies lower than ~ 10 meV, Tic becomes relevant and 
trigonally warps the band structure, which now presents 
four Dirac cones (see Fig. [2]). One of the Dirac points 
(D) remains at the center q = 0, whereas three addi- 
tional cones (A, B, and C) are arranged in a triangle 
around the first one. 



The effective Hamiltonian 

■Ht = Hfa + -He + Ha 



(8) 



was introduced in Ref. |l5| to take into account a trans- 
lational mismatch between the two graphene layers. A 
more microscopic discussion of the model may be found 
inRefs. [iiandlii. 



2. Rotational default 

In the case of a twisted (or rotationally-faulted) bi- 
laycr, lattice-inversion symmetry is broken. For small 
and moderate twist angles, the model Ht + Ha is an 
approximation that yields the correct shape of the en- 
ergy spectrum and the right topological properties of 
the original system, such as the degeneracy of the zero- 
energy Landau levelJ^ The full band structure requires 
taking into account the commensurability between the 
rotated layers and the resulting Moire patterns.— Trig- 
onal warping within the twisted bilayer system is likely 
to be negligible since the orbital mismatch renders all 
hopping parameters small compared to t or t±, such that 
Tic = 0. Notice furthermore that also t± is significantly 
lowered by the twist. For this particular reason, we do 
not consider Ht in Eq. ([S]) as a universal Hamiltonian 
for bilayer graphene but rather as a model that correctly 
interpolates between several configurations that exist un- 
der various experimental conditions. 

Notice that interaction effects generate the same dis- 
tortion Ha, both witbii and without^ trigonal warping. 



C. Band Structure 

The band structure of Ht in Eq. ^ is plotted in Fig. 
[3] for various values of A. 



1. Slide and strain deformation 

The translational and strain constraints may be ac- 
counted for by adding a constant shift 



Ha 



-A 
-A* 



(6) 



to Hamiltonian ([5]) that represents the only relevant per- 
turbation whenever time-reversal and lattice-inversion 
symmetries are preserved)^ Hence translation and strain 
constraints inevitably give rise to the term ^. For in- 
stance, a small sliding deformation renders the t' hop- 
ping anisotropic due to different orbital overlaps. In 
a similar fashion to the anisotropic honeycomb lattice 
problem )^^i^^ the renormalized amplitude modifies the 
continuum approximation by shifting the momentum by 
a constant value A, 



CTT — > CTT ■ 



A. 



(7) 



Undistorted case 



Without any distortion (A = 0, see Fig. [2), the band 
structure is trigonally symmetric, with a central cone D 
and three peripheral ones, A B and C positioned at 



^ = (0,0), A 



B/C^\-,,± 




26' 2& 



(9) 



within a valley. A Taylor expansion of the energy disper- 
sion around the four Dirac points yields 



Eoiq) = cJql + q: 



EA{ci)^cJql + %l, 



Eb/c{ol) = cJlql + 3g2 ± ^Viq.qy, (10) 
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FIG. 2. (Color online) (a) Band structure of the perfectly 
AB-stacked bilayer graphene around one of the valleys, (b) 
Position of the remarkable points in reciprocal space. We 
label the four Dirac cones from Ato D and the corresponding 
saddle points AD, BD and CD. While the Dirac points all 
reside at zero energy, the saddle points have an energy E — 
c? /Ah. The wave vectors are measured in units of c/h and the 
energy in units of <? /b. 



such that one may define averaged Fermi velocities, that 
is vd = JVxVy = c for the D cone and va ~ vb ^ vc = 

-\/3c for tire satellite ones4 Three saddle points join each 
peripheral cone to the central one, see Fig. [21 and arc 
located at 



^"={-21.'" 



™/5B^U.±f 



fll) 



As a consequence of the trigonal symmetry, they occur 
all at the same energy 



^AD ~Es ~ — 



c 
46' 



that is 



^BD ~^CD ^^s' - ^, 



Es^Es' = \(j] tx=^lmcV, 



(12) 
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FIG. 3. (Color online) Band structure of bilayer graphene 
around the K valley for the Hamiltonian "Ha, with a real 
value of A. For A < 0, the two cones D and A start to 
merge [panels (a), for A = —0.1c /b], and give rise to a local 
minimum after the merging transition, at A = —c?/Ab, [panel 
(b) for a value of A = —0.32c /b\. The opposite case of A > 
[panels (c) for A = 0.3c^/b and (d) for A = 0.92c^/b], reveals 
the merging of three cones at a time, D B and C, or triple 
merging. The wave vectors are measured in units of c/b and 
the energy in units of c? /b. 



Dirac cones are then moved from the positions ^ to 
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26' V462 6 ' 



(14) 



where we have used the values of Eq. ^ for the hopping 
amplitudes. 



with the averaged Fermi velocities 



2. Deformation along an axis of high symmetry 

The trigonal point-symmetry is broken as soon as A 7^ 
0. We first consider the case of a real-valued constant, 
corresponding to an applied deformation along the y-axis 
[see Fig. [He)], that is an axis of high symmetry. The 



Vc2 + 46A (2c - Vc2+46A ) 



v\^ V<?+46A(2c+ Vc2 



46A 



n|/P = v/(5c2 _ 45^)2 _ i6c2(c2 - 6A). (15) 



Moreover, the positions of the saddle points arc shifted 



from those described in Eq. (|TT|) to 



AD^[--.0 



(16) 



BD/CD 



3c2 + 46A 
126c ' 



±^ V(3c2 _ 4&A)(9c2 + 4feA) ) 
12oc y 



The saddle points BD and CI? are at the same energy 
1 



Es' = 



12V3&( 



■(3c2-46A)3/2, 



whereas that between AD is found at 

Ss-- + A^i;s,. 

46 



(17) 



(18) 




FIG. 4. (Color online) In the plane [E, A), the positions of 
the two saddle points Es and Egi define four distinct regions. 
The sector L denotes energies that are below both saddle- 
point energies Es and Es' , whereas H describes energies E > 
Es,Es'- The sector M is defined as energies E, with Es < 
E < Es', and M' for Es' < E < Es. Here the energies are 
given in units of (? /Ab. 

Whenever Es ^ Es', it is possible to define different 
low-energy regions that make a distinction between the 
four cones (Fig. S]). These regions turn out to be useful 
for the discussion of the LL spectrum in the Sec. IIIII 
For instance, in the case A < 0, Es < Es' and the A 
and D cones start to move closer together and eventu- 
ally merge when Es = 0, that is at A = — c^/46. The 
other two cones stay apart, well separated [see Fig. ^a.)]. 
Exactly at the merging transition, the band dispersion is 
a semi-Dirac one, quadratic in one direction, linear in the 
other, whereas beyond the transition a gap opens with a 



quadratic dispersion in both directions [Fig. |3ljb)]. We 
emphasize that this merging transition between the pair 
of Dirac cones is exactly the same as in the case where the 
two Dirac cones were related by time-reversal symmetry, 
as discussed in Refs. [22| and |23|. 

On the other hand, when A > one has Es' < Es, 
such that the D, B and C cones converge to a common 
point whereas A stands alone. The three cones are cou- 
pled at energies around Es' [see Fig. |3l^c)]. At the (triple) 
merging transition, A = 3c^/46, the crossing bands bear 
a complex boomerang shape [Fig. El^d)], while further in- 
crease of A does not open a gap in the band structure, in 
contrast to the above-mentioned merging transition for 
A < 0. This difference may be understood in terms of 
winding numbers that play the role of topological charges 
and that are described in detail in Sec. IIIDI 

Notice that, since trigonal warping is a structure at 
very low energy (^ 10 meV), a small perturbation is 
sufficient to drive the system into one of the merging 
scenarios. For instance, a (triple) merging of the Dirac 
points occurs at a very small (^ O.IOA) displacement of 
one graphene layer with respect to the other one, where 
we use perfect AB stacking as the referenceJ^ 



3. Deformation along an unspecified axis 

In addition to the above distortion along a high- 
symmetry axis of the lattice, we consider the more gen- 
eral deformation along an arbitrary axis which corre- 
sponds to a complex- valued A. Fig. \5\ shows the evo- 
lution of the position of the Dirac points when increasing 
A for different values of the angle 9 defined as A = |A|e*^. 
The complex position of the merging point in reciprocal 
space is 



.(0) = ^e^''-^^) 



(19) 



where the angular dependence of the merging angle 
i9,„(6') is given by 



tan0 



2 sin -dm ~ sin 2iD„ 
cos 2'dm + 2 cos "dn 



(20) 



and is plotted in Fig. IH^a). For a given angle 9, the 
merging is reached for a critical value Am{9) given by 



A™(0) = ^[5 + 4cos37?„(0)]l/^ 



46 



(21) 



which is shown in Fig. IHlJb). 

Eq. (|20l) and Fig. IDJa) reveal that the angular depen- 
dence of the merging point, i?„i, is not linear in the angle 
of the deformation axis, 9. This is best captured in the 
vicinity of 6' = or 27r/3 where the slope of the Fig. [D^a) 
increases abruptly. Most saliently, a slight deviation from 
a high-symmetry axis (9 = 0, ±27r/3), i.e. an infinites- 
imal imaginary contribution to the shift A, renders the 
triple-merging point unstable. Indeed, as one may see 
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FIG. 5. (Color online) Motion of the Dirac points in mo- 
mentum space for different values of the displacement angle 
e. The units are such that c = 1, & = 0.06. When A = 0, 
the central Dirac point is surrounded by three Dirac points at 
distance c/b (red dots). When increasing A, the central Dirac 
point merges with one of the three Dirac points, leaving the 
two remaining points isolated. When varying d, the position 
of the merging point draws a circle of radius c/26. The full 
curves represent the positions of the Dirac points until two of 
them merge. The dashed curves represent the position of the 
remaining Dirac points after merging of the other two. 



from Fig. [5l only two Dirac points merge, whereas the 
third one remains isolated. From this perspective, one 
can qualify the triple-merging scenario as unstable. How- 
ever, one can argue that for moderate residual chemical 
doping or in the presence of disorder, the difference be- 
tween a triple-merging points and a single-merging point 
with a close-by extra Dirac cone is smeared out, such that 
the study of the triple-merging scenario may still provide 
physical insight. 
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FIG. 6. (Color online) (a) Dependence of the merging position 
angle i^m as a function of the angle of deformation 6. The 
vertical dashed line for 6 = 27r/3 indicates the 27r/3 rotational 
symmetry, (b) Polar plot of the angular dependence of the 
critical value Am(S) of the deformation at the merging. A,n 
is given in units of c? /Ah. 



ing or those encountered in twisted bilayer graphene^ 
that are not accompanied by a gap opening. The na- 
ture of the different merging transitions turns out to be 
determined by the underlying topological properties of 
the band Hamiltonian. In this section we discuss these 
merging transitions in terms of winding numbers that 
play the role of topological charges the sum of which 
is conserved across the transitions. Furthermore, these 
winding numbers play an eminent role also in the pres- 
ence of a magnetic field, where they determine the num- 
ber of zero-energy modes and where they intervene in the 
semi-classical treatment that describes to great accuracy 
the LL spectrum obtained from the full solution of the 
quantum- mechanical equations (Sec. lIVp . 

In the vicinity of band-contact points, the system may 
be described in terms of the effective two-band Hamilto- 



D. Winding numbers 



In Sec. IlICi we have encountered different merging 
types of Dirac points that fall into two classes: whereas 
the merging transition is associated with the opening of 
a band gap, there are transitions, such as triple merg- 



H(q) 





h^{q) +ihy{q) 







(22) 



diagonalization of which yields the energy spectrum 



fA(q) = ^\/hl{q) + ^H(q) and the eigenstates 



"^ " V2 V ^e*^- 



(23) 



where tan0q = hy{q)/hx{q), and A = ±1 denotes the 
band index. The relative phase (/)q exhibits a particu- 
lar topological structure that we discuss in terms of the 
pseudospin map, which is defined as 



h:q — > {h^{q),hy{q)}. 



(24) 



Because of the single- valuedness of the wave functions 
([^5]) , the map h — [h^ (q) , hy (q)] must retrieve its original 
value, modulo 27r on a closed path that starts and termi- 
nates on a precise value qo . All closed paths therefore fall 
into distinct homotopy classes that are described by the 
integer wc, which is an element of the homotopy group 
Tri{S^) associated with the map h from the closed path C 
(with the topology of a circle S^) in reciprocal space to 
closed paths in pseudospin space. In order to calculate 
this integer, which is the pseudospin winding number, one 
needs to integrate the Berry connection Aq — i^/'^Vq^ 
over the closed path, in terms of the wave functions ([^3)) 
and the reciprocal-space gradient Vq = {d/dq^^d/dq^). 
One obtains 



"P'-^^i 



qV^q 



dq. 



(25) 



As such, w{C) is nothing other than the Berry 
phase^i^ within a factor tt calculated over the path C. 
However, we avoid the name "Berry phase" in the present 
context for two reasons. First, the Berry phase does not 
necessarily need to be an integer, as it has been shown e.g. 
in the case of gapped graphene (or boron-nitride) where 
the Berry phase explicitly depends on the energy of the 
path)^ Only the topological part, which should then be 
viewed as the winding number, of this Berry phase de- 
termines the chiral properties, such as those revealed by 
the LL spectrum in the semi-classical approach discussed 
below. Second, a quantum-mechanical phase is defined 
modulo 27r, and one would therefore not expect different 
physical properties for ttwc as compared to for even 
values of wc or w for odd valuesi^i^ However, relevant 
properties of the level spectrum, such as the dispersion 
relation at intermediate energie o^'^^ , the degeneracy of 
the zero-energy modes, and their protection, depend sen- 
sitively on the precise value of wq ■ Notice that wq is an 
additive quantity - if one dcvides the surface E enclosed 
by the path C into distinct pieces, E1...S7V, the winding 
number is the sum of the partial ones calculated over 
paths Cj encircling the surfaces S j , 



N 



'iC) = J2wiC,). 



(26) 



Furthermore, if one of the merging transitions discussed 
in Sec. IHCI takes place inside a path C, the winding num- 
ber is a conserved quantity. It is simply the sum of the 



winding numbers calculated on paths around the orig- 
inal band-contact points before the merging transition 
and may thus also be viewed as a topological charge. 

As an example, wc plot the map (|24p in Fig. [7] for the 
Hamiltonian "Ht + He + ^A hi different configurations 
corresponding to the sectors L, H, M and M' of Fig. 01 In 
the low-energy sector (L), for energies below both saddle 
points Es and Es', all Dirac points are resolved, and 
there exist thus closed loops encircling each of the points 
[Fig. El^a)]. The vicinity of the points A, B, and C is 
then described by a charge -1-1 each, whereas the central 
point D carries a charge —1. At energies larger than Es 
and Es' [sector H, Fig. [D^b)] , aU closed loops necessarily 
enclose all points, and the topological charge is therefore 
the sum (-1-2) of all individual Dirac points resolved at 
low energies. This situation is to be contrasted to the 
sector M, for energies E with Es < E < Es', [Fig. [Tie)]. 
The points A and D are then necessarily enclosed by all 
corresponding loops, such that the charge is 0, whereas a 
second class of loops can still resolve the points B and C 
(charge -f 1 each). In the sector M', for Es' < E < Es, 
the three points B, C, and D can no longer be resolved 
(loops of charge -fl), whereas A remains a Dirac point 
with charge +1 [Fig. Uld)]. 




,'t\\\\ -, - -^ 

M ^ It ii \ V \y vvwv 
; T H ^ M \\\yyxv\y 

MUM >i\vy\\\\y\y 

H H It "t 



M H itu H xvy 
} H n H "i >i >, V y 




(c) 






"'ZZ\\\ 



it at*, ma rfi 




M M ii >i >t «t >i "i It M y 
M ^ M "i >i >t K \ V V >i W, 

n It It \ ^xyvx^^vvv 
It <! «t V »l^^^^\xy^yvx^ 

"t >t "ivyy 

•t "i X vw 




i i U * i i * f(l! I'I'I"''""' 




*iii*JJ.i 




WW . . . 
H \ XXX\\ 

H It v "t X xyy 

!!!S{!!555 



FIG. 7. (Color online) Winding numbers for the pseudospin 
map. The square indicates a charge —1 and the circles indi- 
cate a charge -1-1. (a) In the sector L, all Dirac points are 
resolved and described by individual topological charges, (b) 
Sector H, the possible closed loops enclose all Dirac points, 
and the topological charge is thus 2. (c) Sector M, whereas 
the Dirac points B and C are resolved (charge 1), the points 
A and D have merged, such that closed loops yield a charge 
0. (d) Sector M', the triple merging envolves the points B, C, 
and D and closed loops yield a charge 1, in addition to the 
charge 1 stemming from the isolated Dirac point A. 

In view of the different merging transitions, we have 
already mentioned that the pseudospin winding num- 



ber is conserved during such transitions. For merging 
(A < 0) of two Dirac points described by winding num- 
bers of opposite sign (e.g. wjj = —1 and wa = +1) ~ this 
is necessarily the case for Dirac points that are related 
by time-reversal symmetry - the topological charges are 
thus annihilated across the transition, such that the zero- 
energy states are no longer topologically protected. One 
therefore observes the opening of a local band gap that 
is associated with the merging of Dirac points with op- 
posite winding numbers [see Figs, ^a) and (b)]. In the 
case of a triple merging ( A > 0) , the sum of the winding 
numbers is wb + wc + wd = +1, such that any path 
enclosing the point where the Dirac points B, C, and D 
have merged carries a winding number +1 also after the 
transition. The (zero-energy) band-contact is therefore 
preserved and the opening of a band gap topologically 
prohibited [see Figs. |3Kc) and (d)]. 



III. LANDAU LEVEL SPECTRUM 




FIG. 8. (Color online) LL spectrum of the Hamiltonian (|29p 
as a function of 5. The parameter /3 is fixed to 0.06. The 
spectrum is cut into four parts delimited by the thick (green) 
lines which depict the energy of the saddle points in units of 
c/Ib = cyeB. The figure shows the existence of four zero- 
energy modes. Two of them (red) are topologically stable. 
The other two (here only the one of positive energy is marked 
with a dashed curve) acquire a finite energy for sufficiently 
large values of 5. 

The considerations of the previous section on the band 
structure in the absence of a magnetic field yield valu- 
able insight into the LL spectrum, which is formed when 
a perpendicular magnetic field, Be^ = V x A, is ap- 
plied to the graphene layers. In this section, we compare 
the LL spectrum obtained from a numerical solution of 



the full quantum-mechanical problem described by the 
Hamiltonian 'Hb + He + ^A in the presence of a magnetic 
field (Sec. IIII A[) to that calculated within a semiclassical 
approximation (Sec. IIIIB[) . A detailed discussion of the 
LLs in the different energy sectors (L, M, M', and H) is 
postponed to Sec. |Vl 



A. Landau Quantization 

The magnetic field may be taken into account with the 
help of the Peierls substitution (for electrons of charge 

-e) 



H: 



eA, 



(27) 



where A is the vector potential. This allows one to in- 
troduce the harmonic oscillator operators 



a := ^ (n^ - iTly) , 



aT = -|(n,-fzn,), (28) 



with [a,a^] = 1. The magnetic length Ib ~ l/^/eB 
26 nm/ 



' B\T] encodes the size of the cyclotron orbits in 
real space. 

In the presence of a magnetic field, the Hainiltonian 
%b+'Hc + 'Ha reads 



c/Ib 







2/3a2 



2/3at2 - y/2a 



V2at 




(29) 



We have rescaled the energy with respect to c/Ib = 
cy/eB, the characteristic LL energy of the central Dirac 
cone D for (5 = 0, and have also introduced the dimen- 
sionless (J5-field-dependent) shift S ~ A/c\/eB as well as 
the parameter 



bVeB 



(30) 



The quantity /3, which measures the amplitude of the 
trigonal warping in units of the inverse magnetic length, 
is a central parameter in the description of the LL spec- 
trum. It may also be interpreted as the inverse of the 
reciprocal-space distance c/b of a peripheral Dirac cone 
{A, B, C) to the central one {D) and the magnetic length 
Ib- Viewed as an energy scale, it is proportional to the 
ratio between the first excited LL {c\j2eB) and the en- 
ergy Es of the saddle points joining the cones [in the 
absence of a deformation (at (5 = 0)]. From \12\ . we have 



V2eB 

Es 



4^2/3. 



(31) 



Finally, the parameter /3 turns out to describe the role of 
magnetic blurring that is described in Sec. IIVBI Notice 
that one might also have performed the Peierls substitu- 
tion in the original four-band model (jlj, as it has been 
done for the case without trigonal warping.~ However, 



the corrections are weak in the low-energy hmit that we 
are interested in, and the effective two-band model ([^^ 
provides a good description of the LL spectrum. 

The numerically obtained spectrum of Hamiltonian 
(|29)) is plotted in Fig. [8] as a function of S, which corre- 
sponds to varying A and/or B as well as the energy of 
the saddle points in order to sustain the same number 
of LL below Es- Notice that we have only plotted the 
spectrum at positive energy, e„ , those at negative energy 
are obtained from the plotted ones simply by adding a 
minus sign, — e^, as a consequence of the particle- hole 
symmetry respected by Hamiltonian (P^ . Increasing the 
value of /3 will then only scroll the levels up in energy 
and reduce the number of modes within the trigonally 
warped area. For this reason, we focus on an arbitrarily 
low value of /? = 0.06. The spectra for other values of (3 
are discussed in Sec. IIVBI 



B. Semiclassical description 

In order to reproduce the spectrum of Fig. [8] and to 
understand the underlying physical properties, we rely on 
a semiclassical analysis. This theory states, according to 



must fulfill 




FIG. 9. (Color online) Semiclassical reconstruction (dashed 
lines) of the spectrum in Fig. [8l for /? — 0.06. The different 
regions of the spectrum are discussed in detail in Sec. [V] The 
blue dots indicate the LL spectrum at the merging (left) and 
the triple- merging transition (right), in which case the LLs 
scale as {n + 1/2)^' '^ and n^'*, respectively, see Sees. IV Gl and 
IV Dl The green lines indicate the energies of the saddle points 
es and eg/ . 



Ac{en) 



k(£<£„) 



(fk==2TTeB{n + j), (32) 



where the mismatch factor 



^=2 



IB 



(33) 



has a contribution 1/2 from the usual Maslov index for 
the harmonic oscillator and a second one, 75, that was 
first identified with the Berry phase^i^ acquired on the 
path C, whereas it has been shown afterwards that only 
the topological part of the Berry phase enters into the 
expressioui^ Here, we express the quantity 75 in terms 



of the pseudospin winding number, 73 
that Eq. dSSl) becomes^! 



wc 



7 



\wc\ 



/2, such 



(34) 



At first sight, the large-n limit of the semiclassical ap- 
proximation could be described in terms of a quantum- 
mechanical Berry phase 7r|wc| modulo 27r, that is one 
identifies all odd and all even winding numbers, if one 
redefines the integer n. However, Eq. p4p bears infor- 
mation about the presence and the number of zero-energy 
modes [see Sec. IIVJ . 

In order to obtain the semiclassical LL spectrum e„ = 
e(n), we numerically invert Eq. ([32]). The results are 
shown in Fig. [5] (dashed lines) in comparison with the 
ones (full lines) obtained from a numerical solution of 
the quantum-mechanical eigenvalue equation [Hamilto- 
nian (pg)) ]. 
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FIG. 10. (Golor online) The different winding numbers at- 
tached to the different pockets imply different quantization 
rules. 



Onsager's argumentr^i^ that the reciprocal-space area 
-4c(en) enclosed by the band contour C, for energy e„. 



In addition to the LL spectrum. Figs. [5] and [5] depict 
the energy of the saddle points rescalcd by the energy 
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c/Ib (thick green lines), 

es = Eslc\feE 



es' 



(35) 



Whenever a LL crosses one of the three saddle points, its 
properties, such as its degeneracy, are drastically modi- 
fied due to the Lifshitz transition involved. Indeed, Figs, 
m and [3] reveal four distinct sectors L, H, M, and M' in- 
troduced in Fig. |4] according to whether the energy is 
below both saddle points [low-energy sector (L)], below 
only one of the saddle points [merging (M) and triple- 
merging (M') sectors], or above both of them for the 
high-energy sector (iJ) [see Fig. [T4la)]. For these distinct 
regions, Onsager's quantization reads differently, because 
the winding number wc takes different values for differ- 
ents types of orbits. The different sectors are shown on 
Fig. [TOland discussed in detail in Sec. El 



IV. ZERO MODES AND SEMICLASSICAL 
QUANTIZATION RULE 

In principle, the semiclassical description is valid at 
large energies, whereas the zero modes require a specific 
(quantum-mechanical) treatment. However, the semi- 
classical analysis and the intervening winding numbers 
provide valuable insight into the degeneracy of the zero 
modes. The discussion of this relation is the issue of the 
present section. In Sec. IIV Al we provide a simplified 
model to illustrate this relation, whereas we discuss the 
LL degeneracy lifting due to magnetic blurring in Sec. 
IVBl 



Relation betw^een zero modes and ^vinding 
number 



We provide here a heuristic argument relating the total 
number of topologically protected zero-energy modes to 
the semiclassical quantization rule. Consider first the 
model Hamiltonian^li^d^ 



H = A 



TftP 

ttP 



(36) 



describing a band contact point with energy spectrum 
e = ±A|q|P and with a winding number p. The global 
parameter A has the physical dimension of an energy 
times the p-th power of a length. In a magnetic field, 
performing the Peierls substitution (|27p and the replace- 
ment in terms of ladder operators (j28p . one obtains the 
LL spectrum in a magnetic field. 



,(B) = ±A(2eB)P/ V"(" - 1) • • • (n - P + 1), (37) 



which, in the large-n limit, may be approximated as 

p/2 



e„(B) ~ ±A 



2eB ( n H 

2 2 



(38) 



This corresponds precisely to the semiclassical quanti- 
zation rule p2p if we identify -p with the total winding 
number wq in Eq. ()34|) . From Eq. p7|) . one notices 
that the p quantum numbers n = 0, ...,p — 1 correspond 
to states at zero energy. Indeed, these states may be 
obtained from the eigenvalue equation 



U 



0, 



which is satisfied for the states 



\n) 



t 



(0) 



for 



0,...,p-l 



(39) 



(40) 



in terms of the eigenstates \n) of the number operator 
a^a, a^a|?i) — n\n). Because of the orthogonality of the 
states, the zero-energy manifold is p-fold degenerate, i.e. 
the degeneracy corresponds to the total winding number 
Wc = p, as stated above. 

The situation is different when there are more band- 
contact points, with different (local) winding numbers. 
Consider a Hamiltonian describing p massless Dirac 
points with winding number -|-1, situated at the com- 
plex positions a^ in reciprocal space, and p' massless 
Dirac points with winding number —1, at the positions 
/3j . (Notice that band contact points with larger winding 
numbers may be obtained by making several positions ai 
or f3j coincide.) The Hamiltonian can be written as 



-1/ _ A I •'1 



(41) 



with /q = nj=i('''^ ~ Pj)IVi=ii^ ~ ^i) f-iid A is a global 
constant of the dimension energy times the (p + p')-th 
power of a length. The total number of Dirac points, 
and thus, after the Peierls substitution (P7)) . the maximal 
number of zero-energy LL, is Wt = p + p' . However, this 
(p -|- p')-io\d degeneracy of the zero modes may be par- 
tially lifted upon merging of two or more Dirac points. 
In order to find the total number of topologically pro- 
tected zero-energy levels, we thus continuously modify 
the parameters 



and 



Pj ^ 0, 



(42) 



so that /q becomes tt^ 



In a magnetic field, assum- 



ing for example that p > p' , this term is of the form 
\/2eB *{a^ayaP~P , and the associated LL spectrum 
reads 



en{B) = ±A(2eB)""/2nP' ^"(^ - 1) ■ ■ ■ {n - p + p' + 1). 

(43) 
The same arguments as those presented in the discussion 
of the Hamiltonian ((36)) indicate that there are Wp = 
p — p' zero-energy levels that correspond to the quantum 
numbers n = 0, ...,p — p' — 1 . Moreover, the same large- 
n expansion as in the case discussed above yields the 
spectrum 

Wp/2 

, (44) 



e„(S) ~ ±A{2eB)'^nP' 



1 w, 
nH i 

2 2 
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which may be cast into the semiclassical quantization 
rule with a winding number wc = Wp = P — p' ■ Again 
the total winding number is identical to the number of 
zero-energy modes. 

These arguments show that, although the maximum 
number of zero-energy LLs is wt ~ p + p' , & quantum- 
mechanical coupling between them partially lifts the de- 
generacy, but Wp — \p — p'\ zero modes remain topologi- 
cally protected. Applied to the model dS]), one has p ~ 3 
and p' = 1, so that the maximal number of zero modes is 
Wt = 4 and the number of topologically protected modes 
is Wp = 2. 



B. Magnetic blurring 

It is apparent from Fig. [9] that the semiclassical treat- 
ment p2p based on Onsager's quantization rule provides 
a reliable description of the LL spectrum in the major 
part of the parameter range. However, it is challenged in 
the vicinity of the saddle points £5 and es' ■ Intuitively, 
one may understand the failure of semiclassical quanti- 
zation if one considers the topological winding number 
pS]) that is calculated from closed loops around the re- 
markable points in reciprocal space [see Fig. [7] - in the 
presence of a strong magnetic field, these loops are at 
odds with quantum mechanics because the components 
of the wave vector are no longer good quantum num- 
bers, such that the images of the loops defined by the 
maps (|24p are constrained by a Heisenberg uncertainty 
relation. 

Indeed, in the presence of a magnetic field, the mo- 
menta no longer verify the simple commutation relation 
[n^cllj,] = but rather a Heisenberg algebra 






(45) 



An immediate consequence of this non-commutative ge- 
ometry is that reciprocal space is now "patched" or 
"blurred" by irreducible regions of area 1/1% below which 
it is impossible to resolve the physical properties of elec- 
trons in a magnetic field. This is similar to the phase 
space of a one-dimensional quantum-mechanical parti- 
cle, which is devided into minimal regions of the size of 
the Planck constant h = 2tt, below which the physical 
properties of the particle cannot be resolved. As a con- 
sequence, the winding of the pseudospin vector cannot 
be determined by paths the area of which encloses less 
than the minimal area of ex 1//^, which plays the role 
of the Planck constant in reciprocal space. If we were 
to define a winding number in the non-commutative re- 
ciprocal space, we should then consider larger and larger 
contours as the magnetic field increases since 1/1% ex B, 
as shown in Fig. [TTJ 

One may thus pictorially understand that whenever 
the irreducible area of 1/1% becomes too large, the rele- 
vant winding contours enclose inevitably more than one 
singularity. This is shown in Fig. [TT] where at low fields 
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FIG. 11. (Color online) Magnetic blurring for zero modes. 
Each contour encloses a minimal surface of ~ l/'s ~ eB in 
reciprocal space (red areas), (a) At low magnetic fields, the 
blurring is low, and each minimal surface contains a single 
Dirac point. The semiclassical quantization rule holds for 
each contour encircling a Dirac point, (b) When the field 
increases, the contour associated with the zero mode encloses 
a larger surface so that the individual Dirac points are no 
longer resolved. The effective winding number experienced 
by the electron is u;p = ^^ w^. 



the contours Ci around the individual singularities en- 
close each a winding number w^ = ±1 whereas at high 
fields the blurred contour encloses a winding number 
^^ Wi ~ 2. Therefore, for a sufficiently strong magnetic 
field, the only relevant quantity is the total winding num- 
ber around all the singularities. 



E^ 



(46) 



as opposed to the total sum of the winding numbers 

Wt=Y,\w,\. (47) 

i 

which sets the total number (but not necessarily pro- 
tected) of Dirac points. In the model Hamiltonian dis- 
cussed in Sec. IIV Al this magnetic blurring may be 
viewed alternatively as an effective merging of the band- 
contact points (|42|) . 

As a consequence of the above arguments, increasing 
the magnetic field induces, even at zero energy, a Lifshitz 
transition that is characterized by a partial degeneracy 
lifting of the zero-energy LL from wt (per spin and valley) 
to Wp < wt, while Wt — Wp levels disperse as a function 
of B. 

In the case of the Hamiltonian (P^ , we generally have 
four Dirac points in the vicinity of (5 = such that one 
expects a four-fold degeneracy of the zero-energy mode 
for low magnetic fields, in agreement with Eq. (|46p . Be- 
cause of the parameter /3 a VB [see Eq. (pO)) ]. the low- 
field limit corresponds to small values of /?, such as in the 
case of the value 0.06 chosen to calculate the spectrum in 
Figs. |5]and|ni Indeed the fourfold degeneracy of the zero 
mode is lifted only when approaching the merging tran- 
sitions from 4 to 2, where /3 diverges as a consequence of 
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that read, to lowest order in Ig, 
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FIG. 12. (Color online) LL spectrum of the Hamiltonian (|29p 
as a function of 5. The parameter jS is fixed to 0.2 (a) and 
0.45 (b). (c) LLs as a function of the parameter /3 in the 
absence of distortion (5 = 0). It exhibits the lift of the fourfold 
degeneracy when the energy of the LLs becomes larger than 
the saddle point energy Es (blue levels). Two zero energy 
levels stay stable while two other levels get a finite energy 
when/3 > 1/4^2^0.18. 



the decreasing reciprocal-space distance between some of 
the Dirac points. However, this degeneracy can also be 
lifted exactly at (5 = by increasing the value of /3, where 
Fig. [TTJand Eq. (|Tf)) indicate that the degeneracy of the 
zero-energy mode is 2 above a certain magnetic field. The 
quantum-mechanical LL spectra, for larger values of /3, 
are depicted in Fig. [T2la) for /3 = 0.2 and for 0.45 in 
Fig. [T^b). as a function of 5. In both cases, one notices 
that the fourfold degeneracy is indeed lifted for all values 
of 5, in agreement with the expectation from magnetic 
blurring. The effect is also apparent in Fig. [T2lc) , where 
we have plotted the (5 = LL spectrum in units of the 
saddle-point energy Es as a function of /3. Indeed two 
branches of the small-/3 zero-energy mode float away ~ 
due to particle-hole symmetry, one increases in energy 
while the other one decreases - while two other branches 
are topologically protected and remain at zero energy. 

Notice that the magnetic blurring in reciprocal space 
may also be understood as a blurring in energy. Indeed, 
the commutation relations (|^5|) induce, via the maps (P^ . 
commutation relations for the pscudospin components 



[hx^hy] 



dhx dhy 

1% [dUydUx 



dhx dhy 

du^ an, 



(48) 



In the vicinity of a Dirac point j with linear band dis- 
persion and a characteristic (possibly anisotropic) Fermi 
velocity {vx.j,Vy,j), the commutation relations (|48p thus 
induce a Heiscnberg uncertainty relation AhxAhy ^ 
Vx,jVy,j/lB = ^j/^B that is precisely on the order of the 
energy gap between the zero-energy level and the first ex- 
cited one. In this picture, the topological winding num- 
bers and thus the level degeneracies associated with indi- 
vidual contours around the Dirac points arc well-defined 
as long as the energy uncertainty y/Ah^Ah^ ^ i^j/^B is 
smaller than the saddle point Es. This argument agrees 
with the expectation that the zero-mode degeneracy is 
lifted once 



Es < V2c\ 



^ 13 



> 



1 



4^2' 




^ '' ^ ^ i' i^W^^^ 






(49) 
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FIG. 13. (Color online) Magnetic blurring for contours at 
higher energy, (a) At low magnetic fields, the semiclassical 
quantization rule holds for each contour encircling a Dirac 
point, (b) When the field increases, the energy contours 
become blurred and tunneling to trajectories enclosing two 
singularities becomes possible in the vicinity of the saddle 
points. The effective number experienced by the electron is 
W(c+C') = wc + wc = for the upper bound of the contour, 
whereas it is |toc| + l^cl = 2 for the lower bound. 

In addition to the zero-energy modes, magnetic blur- 
ring also plays a role in the degeneracy lifting of higher- 
energy LLs in the vicinity of the saddle points, where the 
semiclassical approximation does not accurately describe 
the LL spectrum [see Fig. |9]. Indeed, the degeneracy 
lifting in the semiclassical approximation is abrupt be- 
cause of the abrupt change in the winding number: for 
energies just below the saddle points, one has discon- 
nected energy contours C and C that become connected 
by a contour C + C for infinitesimal energies above the 
saddle points. However, this abrupt transition is blurred 
because not only the above-mentioned smallest contours, 
which are responsible for zero-energy modes, need to en- 
close a minimal surface of ~ l/'s' but also two contours 
corresponding to successive energy levels (see Fig. [T51) . 
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The resulting uncertainty about whether a contour in 
the red region in Fig. [T3] is connected or disconnected 
yields an uncertainty in the winding number, such that 
the variation of the LLs in the vicinity of saddle points 
is smoother than that expected from the semiclassical 
analysis. 



V. DETAILED ANALYSIS OF THE SPECTRUM 

The semiclassical and topological theories presented in 
Sees. IIIIBI and IIVI respectively, allow us to discuss in 
detail the different properties of the LL spectrum in Fig. 
[51 From the semiclassical quantization p2p with appro- 
priate values of the winding number wc^ we obtain the 
semiclassical spectrum in the different energy regions sep- 
arated by the saddle-point energies. Unless stated explic- 
itly, we discuss only the orbital degeneracy for a single 
valley {K or K') and a single spin - the full degeneracy 
is then given by the orbital degeneracy times the fourfold 
spin- valley degeneracy. 



A. Undistorted Case 5 — Q 

We begin our reconstruction of the spectrum by plot- 
ting the LL at (5 = (A = 0), that is without any distor- 
tion, see vertical lines in Figs. [Hfa) and (b). 

As long as the energy satisfies e < eg = eg' j the LLs lie 
in the central region of the low-energy sector (L), see Fig. 
[TW a). In the absence of a magnetic field, the low-energy 
spectrum is that of Fig. ^a) and consists of four Dirac 
cones, A to D. The four Dirac cones give rise to four 
disconnected Fermi pockets the area of which is 



Aq 



TT^ Ap 



3c2' 



(50) 



where Ao is the area of the central cone (D) and Ap 
that of the other peripheral ones. In the vicinity of 
each isolated Dirac cone, one has a topological charge 
of wc = ±1, as in the case of monolayer graphene, such 
that Onsager's quantization rule (|32p yields 



27r n 



1 



\wc\ 



eB = 2TmeB = A{E) (51) 



=> Eo{n) = cV2neB Ep{n) = cVGneB. 

At zero energy and for moderate magnetic fields (/3 = 
h^feBjc <^ 1, as in Fig. [5]), one thus obtains a fourfold 
orbital degeneracy of the zero-energy level because each 
of the four Dirac cones yields an n = LL. Therefore 
the total degeneracy of the zero-energy level, taking into 
account again the additional fourfold spin- valley degen- 
eracy, is 16. From the topological point of view, this 
is related to the presence of four well-separated Dirac 
points that are not yet blurred by the relatively moder- 
ate magnetic field, /? ^ 1, in such a manner that the four 




FIG. 14. (Color online) Different sectors of the LL spectrum, 
for /3 = 0.06. (a) Landau levels in the low-energy sector (L) 
obtained from the semiclassical approximation A = 2-neBn. 
The red levels correspond to the quantization of the D cone, 
the blue dashed levels to the quantization of the A cone. 
The last set (purple dashed-dotted) of levels is twofold de- 
generate since it corresponds to the quantization of the two 
cones labelled B and C. (b) Landau levels in the high- 
energy sector (above the saddle points) . The blue levels in the 
M region are obtained from the semiclassical approximation 
A = 2neB{n + 1/2) corresponding to the absence of wind- 
ing number. The purple levels in the M' zone are obtained 
from the condition A = 2-neBn of the pocket issued form the 
B, C, D cones. The blue dots indicate the LL spectrum at 
the merging (left) and the triple-merging transition (right), 
in which case the LLs scale as (n -I- 1/2)^'^ and n'''^, respec- 
tively. In the high energy region H, the four Dirac pockets 
have merged into a single pocket with a total winding num- 
ber 2, so that the red levels are obtained from the condition 
A = 2TTeB{n+l/2). 



related winding numbers are decoupled and give rise to 
four zero-energy LLs. 

For higher (relativistic) LL, we find three times more 
LLs for the central cone because Aq = 3Ap. This ex- 
plains the fourfold degeneracy displayed one time out 
of three in Fig. [Tif aV Indeed, every third LL is not 
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only associated with the central Dirac point but also 
with the three peripheral ones; hence its (accidental) 
fourfold orbital degeneracy, whereas the other LLs arc 
non-degenerate because they are associated with the cen- 
tral one only. As we have previously discussed [see Eq. 
(fTO|) ]. the averaged Fermi velocity at the points A, B, 
and C is v3 times larger than that, c, around the D 
point, VA = vg/Q = v3v£). The three times denser LL 
spectrum associated with the central cone is therefore a 
consequence of this relation between the Fermi velocities 
and of the approximate LL dispersion 



Vi I — 



(52) 



around the Dirac points, with their typical ^/n scaling, 
for 2 = A, B, C, or D. 

Crossing the transition lines at e = £5 = es' , the spec- 
trum undergoes a transition to the high-energy sector 
(H). sec Fig. [Tif b). At zero magnetic field, the three pe- 
ripheral Fermi pockets merge with the central one. This 
change in the topology of the Fermi surface has conse- 
quences for the degeneracy, in the sense that Onsager's 
quantization rule indicates that there is only one set of 
LLs associated with the simply connected Fermi surface. 
There is thus no more orbital LL degeneracy and, besides 
the trigonal deformation of the Fermi surface, the band 
structure is approximately parabolic such that the LLs 
scale as n, as in the usual description of Bcrnal-stacked 
bilayer graphene {^Jn{n — 1) ^ n — 1/2). 



B. Slightly distorted case < |5| < 1 

For small non-zero values of the parameter S, the sad- 
dle points occur at two different energies, es 7^ ^S'j a-nd 
one therefore needs to distinguish three different energy 
sectors, that is L, M, and H for (5 < and L, M', and H 
for 5 > [see Fig. [TU]. Below the energies es and es', 
the band structure is that of Fig. EKb) or (e), comprising 
four Dirac cones albeit with no trigonal symmetry due 
to finite distortion of the bilayer. For small values of /3, 
the picture obtained in the discussion of the 6 = case 
remains essentially unaltered at zero energy. The pres- 
ence of four distinguishable Dirac cones yields a fourfold 
zero-energy level that is insensitive to the slight geomet- 
ric deformation of the perfectly trigonally-warped case. 
As a consequence, the zero modes remain untouched over 
a wide range of 6 distortion around 0. 

On the other hand, the higher LL are not topologically 
protected and the breaking of the trigonal symmetry in- 
duces an immediate lift of orbital degeneracy, as pictured 
in Fig. [131 Our choice of a real valued parameter A im- 
plies that the Dirac cones B and C are related by mirror 
symmetry, such that vb = vq- The corresponding LLs 
[thick blue lines in Fig. [Mia)] are therefore twofold de- 
generate and experience the strongest decrease in energy 
with increasing 5 > Q because their average Fermi veloc- 
ity is decreased [see Eq. ([TS]) ]. The other two sublevels 



have a single orbital degeneracy, corresponding to the 
Dirac cones A and D. As one may see from Eq. ([T5|), the 
Fermi velocity of the central cone D decreases (quadrat- 
ically in J), whereas that of the cone A increases linearly 
in 5. As a consequence, the energy of the LLs associated 
with D [red lines in Fig. [TiT a)] is decreased both for pos- 
itive and negative values of (5, whereas the A-conc LLs 
increase linearly in energy with increasing 5. 

Above both saddle points, i.e. in the sector H, varying 
5 always yields a decrease in the size of the unique Fermi 
surface so that the non-degenerate LL are enhanced in 
energy, as one may see in Fig. I14f b). 



C. Merging transition 5 <§; — 1 



In Sec. lIVBi we have shown that the fate of the zero- 
energy level is determined by the parameter /3 ~ upon in- 
crease of /3, one obtains a magnetic-field-induced Lifshitz 
transition from a fourfould degenerate to a twofold de- 
generate level. Whereas this picture is roughly the same 
for small values of |(5|, it needs to be modified when ap- 
proaching the zero-field merging transition, that is when 
the saddle point energy es vanishes. As one may see 
from Fig. [HJ one notices a significant departure from the 
semi-classical approximation. The A and D Fermi pock- 
ets merge indeed into a single one and the corresponding 
orbital degeneracy of the LLs is changed. Indeed, be- 
cause of the decrease in energy of the saddle point £^5, 
the latter is only higher in energy than the typical scale 
\/2c\feB for the separation between the lowest LLs if 



1< 



1 



4V2/3 






(53) 



This is a generalization of the criterion (|49p for the undis- 
torted case (5 = 0. Based on the criterion ([5^ . one there- 
fore expects the zero-mode degeneracy to be partially 
lifted at (5 ~ \f2 — 1/4/3 [that is (5 ^ —3 for our above 
choice /3 = 0.06], in good agreement with the numerical 
results depicted in Fig. [8l Directly at the merging tran- 
sition, that is for (5 = — c/46veB, one obtains a LL spec- 
trum with levels that scale as e„ = 2Afi^^'^(n + 1/2)^/^ 
with A = 7r[3/r(l/4)]2/3 ~ 1.173, in agreement with the 
merging transition of Dirac cones with opposite Berry 
phases 1^1^ Upon a further decrease of i5, the merging of 
the cones A and D is associated with a gap opening [see 
Fig. [3)Jb)] such that the corresponding LL spectrum is 
shifted to higher energies, as may be seen on the left-hand 
side in Fig. [T4l b). Apart from the shift to higher ener- 
gies, these LLs corresponding to the merged points scale 
linearly in the LL index n, as one expects for parabolic 
bands {A^E) oc E (x n). Because of the distance in 
energy from the saddle points, the semiclassical approxi- 
mation agrees well with the numeric spectrum, as can be 
checked in Fig. [Hi 
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D. Triple-merging transition 5^1 

In the opposite limit, for 6 > 0, the A cone remains 
apart and its Fermi velocity va is increased. The en- 
ergy of the LLs is therefore enhanced and well described 
within the semi-classical approximation, as may be seen 
in Fig. [3] The A cone is unaffected by the transition line 
indicating the saddle point eg' because it is not involved 
in the triple-merging process, as opposed to B, C and D, 
that form a boomerang-shaped Fermi surface. The latter 
become coupled through the Lifshitz transition at Es' , as 
can be observed from the departure from the semiclassi- 
cal approximation in Fig. [S] Before this transition, all 
three Fermi pockets increase in size, with a higher rate 
for B and C than for D, such that the twofold degenerate 
LLs corresponding to the points B and C decrease faster 
in energy than those of the central cone D. 

There are only two zero-energy LLs since the total -1-1 
topological charge of the boomerang pocket gives rise to 
a unique topologically protected mode. Equivalently, the 
magnetic field has reached such a value that the total 
winding number Wp is the only relevant quantity. 

Precisely at the triple merging point {Es' = 0), the 
LL scale as n^^^ [blue dots in Fig. |Hl[b)], as far as our 
numerical accuracy is concerned. After the triple merging 
transition, the LL oi D = B = C scale with a different 
power law, almost linear in n. Increasing i5 increases 
the energies of all the LL because of a decrease of the 
combined-orbit area. 



E. LL spectrum for an imaginary value of S 
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FIG. 15. (Color online) An imaginary value of A (S = ""/S) 
induces a merging between the B and D cones [panel (a) for 
Im(A) > 0) or between the C and D cones [panel (b) for 
Im(A) <0]. 

For the sake of completeness, we present here a LL 
spectrum for a purely imaginary A in the case of a de- 
formation in the x-axis {9 — 7r/2 in Fig. [HI sec Sec. Ill CI 
3). Then, the two directions of deformation [±Im(A)] 
are equivalent, as pictured in Figs. WL^) andfTSl so that 
the saddle point energy and the LL spectrum is now sym- 
metric in Im(A), as seen in Fig. [121 When Im(A) > 0, 
the cones B and D merge leaving the cones A and C 



isolated, whereas for Im(A) < 0, the role of the B and 
C points is interchanged. In a magnetic field, one distin- 
guishes the LL sequence from the four cones below the 
saddle point energy, as well as the Lifshitz transition near 
the saddle points. Notice that the accidental degeneracy 
of the B and C cones, which we have encountered for 
real values of A, is now lifted and that the LLs associ- 
ated with the A cone are symmetric in Im(A) (Fig. I16p . 
Indeed, the A cone remains isolated and does not take 
part in the merging transition for any value of Im(A). 
Its LL spectrum therefore remains relativistic with the 
typical y/n scaling. In the merging sector, for energies in 
between the two saddle points, the C cone provides an 
additional set of relativistic LLs for Im(A) > 0, whereas 
this set is provided by the B cone for Im(A) < 0. As 
in the case of the merging transition discussed in Sec. 
IVCl beyond jA^j ^ (c2/46)(6^/3 - 9)^/2 ~ 1.18c746 
{\5,n\ ^ 4.92, for /? = 0.06 as shown in Fig. ^ the 
merged cones [B and D for Im(A) > 1.18c^/46 or C 
and D for Im(A) < — 1.18c^/46] are accompanied by the 
opening of a local gap the ^-dependence of which is indi- 
cated by the thick blue line in Fig. [THJ Consequently the 
associated LLs are non-relativistic with a linear-n scal- 
ing because of the annihilation of the winding numbers 
of the merged cones. 




-5 5 

Im[5]=Im[^\lc^reB 

FIG. 16. (Color online) Landau level spectrum for 9 — ■k/2 
(/3 — 0.06). In this case, the spectrum is symmetric in the 
displacement lm(S) — Im(A)/cveB. The saddle point ener- 
gies are indicated by the green curves and the energy of the 
local gap beyond merging is indicated by the blue curves. 
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Consequences for magneto-transport 
measurements 



We finish this section with a brief discussion of the 
consequences of the above picture for magneto-transport 
measurements, namely in the context of Hall quantiza- 
tion. Such experiments have been performed both in bi- 
layer graphene in the low-energy limit ^ as well as in 
samples with a twist between the two layers 42*^ Remem- 
ber that the model ^ investigated above also accounts 
for the case of twisted bilayer graphene, with moderate 
twist angles, if one sets c = or in the limit c^/6 <C A. 
In this case, A is a function of the twist angle. 

In all measurements, an eightfold degeneracy of 
the zero-energy level, with no additional quantum-Hall 
plateaus in the range —4 < v < A, has been observed. 
This indicates, in addition to the usual fourfold spin- 
valley degeneracy, a twofold degeneracy of the zero- 
energy level. In the case of twisted bilayer graphenej^i^ 
this is an indication for the twofold topological degener- 
acy associated with two Dirac points characterized by a 
unit winding number with the same sign»ii For a bilayer 
sample with no twist fi^ the observed eightfold degener- 
acy indicates a prominent non-zero value of A since one 
would expect, based on the above arguments for A ^ 0, 
a 16-fold degeneracy of the zero-energy level, i.e. no 
quantum-Hall plateaus in between —8 < v < S. It has 
been argued that the rather large value of A cannot be 
explained by strain (or a displacement of the two lay- 
ers) alone and that interaction effects are likely to play 
an important role^^ in which case A plays the role of a 
(nematic) order parameter)^ 

In higher LLs in trigonally-warped bilayer graphene, 
the degeneracy depends both on the value of the saddle 
point energy '^ |A| as compared to the magnetic energy 
scale cVeB, as well as on the phase of A = |A| exp(i6'). 
In the low-energy sector (L), wc have shown that for a 
real value of A (0 = or tt, modulo 27r/3), the Dirac 
cones at B and C are related by mirror symmetry and 
their LLs are thus (2 x 4)-fold degenerate in the low- 
energy (L) and merging (M) sectors. One would there- 
fore expect a jump of Az/ = 8 in the Hall conductance 
whenever the Fermi level crosses such a level, whereas 
the LLs associated with the points A and D are only 
spin-valley degenerate, associated with a jump Az/ = 4. 
Notice that the mirror symmetry is immediately broken 
in the case of a non-zero imaginary part of A, i.e. when 
^ 7^ or TT (modulo 27r/3), such that all LLs arc then 
fourfold spin-valley degenerate only. This fourfold de- 
generacy is also the generic case in the other energy sec- 
tors [M' and H). Experimentally, quantum- Hall features 
have been observed at i^ = ±4, ±8, ±12, ...,— such that 
the LLs are only spin- valley degenerate. Whereas this 
sequence is identical to that of bilayer graphene without 
trigonal warping, the vS-scaling of the gap between the 
zero-energy level and the first excited one indicates that 
the low-energy sector is nevertheless governed by Dirac 
cones with a linear dispersion relation, as one would ex- 



pect in the sectors M and M'. 

The situation is different in twisted bilayer graphene, 
where one expects eightfold-degenerate LLs below the 
saddle point at E ^ |A|. whereas they arc fourfold- 
degenerate at E > |A|. Since the value |A| can be 
tuned to great extent by the twist angle, one may ex- 
pect to see this crossover more easily than in trigonally- 
warped bilayer graphene. From an experimental point 
of view, Lee et al^ investigated an epitaxially grown 
sample on SiC, with typical twist angles of 2.2°. In 
this sample an eightfold degeneracy of the zero-energy 
level has been observed, whereas higher LLs are fourfold 
spin- valley degenerate, that is a filling-factor sequence of 
V = ±4, ±8, ±12, .... Sanchez- Yamagishi et al. have in- 
vestigated a twisted bilayer sample fabricated by PMMA- 
transfer technique of two monolayer samples on hexa- 
boron-nitride4^ In this case, the observed sequence of 
quantum-Hall plateaus is i^ = ±4, ±12, ±20, ..., that is 
eightfold-degenerate Landau levels also at higher energy. 
This indicates a large value of | A | and of relatively large 
twist angles. 



VI. CONCLUSIONS 

In conclusion, we have investigated a continuum model 
that accounts for the presence of two and more Dirac 
points in the dispersion relation. The model describes the 
low-energy physical properties of bilayer graphene with 
a stacking default, either a translational displacement of 
one graphene layer with respect to the other one, as com- 
pared to the perfectly AB-stacked case, or strainJ^i^ 
Furthermore, this model also accounts for a rotational 
stacking default (twist) if one neglects the linear term 
Tic in the low-energy Hamiltonian ([5])Ji Whereas the 
number of Dirac points is determined by the interplay 
between the different microscopic parameters, the total 
winding number of +2 topologically guarantees the pres- 
ence of at least two Dirac points (with winding number 
+1) or a single parabolic band-contact point (with wind- 
ing number +2). 

In the presence of a magnetic field and Landau quan- 
tization, the winding number yields a doubly degenerate 
zero-energy level that is topologically protected. We have 
studied the LL spectrum in the framework of a semiclas- 
sical treatment and find that it describes accurately the 
numerically obtained one in a large parameter range. The 
semiclassical theory allows for a detailed understanding 
of the LL spectrum in the sense that one may associate 
certain levels with particular Dirac points and determine 
the degeneracy of the levels. Furthermore, the degener- 
acy lifting is understood in terms of connections between 
Fermi pockets. 

However, the semiclassical approximation, which is 
based on a quantization of reciprocal-space orbits and 
the topological charge (i.e. the winding number), breaks 
down in the vicinity of saddle points in the (zero-field) 
dispersion relation as well as at zero energy in the high- 
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magnetic- field limit. The physical origin if this break- 
down is the definition of the topological charges in terms 
of closed reciprocal-space orbits, which change abruptly 
at the saddle points when two or more Fermi pockets 
become connected. Indeed, the definition of topological 
charges needs to be revisited in the presence of a mag- 
netic field that quantizes reciprocal space into patches 
of size ~ 1//^ oc B because the components of the 
kinetic-momentum operator no longer commute. This 
effect blurs the reciprocal-space orbits and smoothens 
the abrupt change in the winding number at the saddle 
points. 

Another effect of this magnetic blurring concerns the 
zero-energy states. Because reciprocal-space orbits need 
to enclose minimal surfaces of r^ l/^s, neighboring Dirac 
points at zero energy are no longer resolved individually 
in the high- field limit. This effect is at the origin of a 



magnetic-ficld-induccd Lifshitz transition, where the de- 
generacy of the zero-energy level, which consists of four 
{wt ~ 4) n = LLs (associated with the total sum wt 
of Dirac points), is partially lifted when increasing the 
magnetic field. Eventually, the degeneracy of the zero- 
energy level is then given by the total topological charge, 
that is iCp = I X^i ""^'Ij ill terms of the zero-field winding 
number Wi — ±1 oi a, single Dirac point. 
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